%Generates instances of autoregressive processes with 
%sparsity patterns in the inverse spectrum which correspond to 
%the two step graph generated form the graph defined by the matrix G.
%Outputs an N time ordered samples of a realization in a matrix x of N by n,
%outputs the auto regressive coefficients A = [I -A_1,..,-A_p], O=Omega
%Equivalent model for white noise innovations B = [B0 B_1,...,B_p] where B0*B0 = Omega
%The AR order p, the variates dimension n.
%The windowed estimate of the covariance matrix C and saves it to a file at 'file'

%clear all
%Set the random seed
%rand_stream = RandStream('mt19937ar','Seed',10);
%RandStream.setDefaultStream(rand_stream);

N = 50;
n = 6;
p = 4;
burn = 100;
%p_e_r = 0.2; %Erdoos reinii probability of an edge being included
%%Generate a random graph for the model by the Erdoos Reinii process
%G  = sparse(zeros(n))
%for i=1:n
%	for j = i:n
%		b = rand();
%		if b<p_e_r;
%			G(i,j) = 1;
%			G(j,i) = 1;
%		end
%	end
%end

%Generate a graph by forming a cycle of n vertices, 
%introduce self edges and edges to the nodes one step ahead

G = eye(n) + diag(ones(n-1,1),1) +diag(ones(n-1,1),1);
G(1,n) = 1;
G(n,1) = 1;
G = sparse(G);

stable = false;
%Iterate untill we get a stable model
while ~stable

	%Generate k random sparse matrices with
	%sparsity pattern equal to that of G and
	%store them into a matrix formed by the
	%concatenation of all teh matrices
	B = zeros(n,n*p);
	for i = 1:p
	 B(:,(i-1)*n+1:(i-1)*n+n) = 1/p*sprandn(G);
	end		

	%Generate the precision matrix for the noise term
	O = sprandsym(G+G');
	
	%TODO: How do we keep it psd independently of n, and not so roughly
	O = O + n*speye(n);
	
	%Calculate B0
	B0 = sqrtm(full(O));
	iB0 = inv(B0);	

	%Generate the regression coefficients
	A  = inv(B0)*B;

	%Form the companion matrix 
	Comp = [-A;[eye(n*(p-1)),zeros(n*(p-1),n)]];
	ei= eigs(Comp); 
	me = max(abs(ei)); 
	fprintf('Maximum eigenvalue of the companion matrix %d\n',me)
	if max(abs(ei)) < 1-eps
		stable = true;
		fprintf('Selected stable model with maximum eigenvalue %d\n',me);
	end
end


%Generate the Ar process
Xt = zeros(p*n,1);

fprintf('Burn in \n')
%Burn in process
X = [];
for j = 1:burn
	Xn = -A*Xt+iB0*randn(n,1);
	Xt = circshift(Xt,n);
	Xt(1:n) = Xn;
end
fprintf('Generating realization \n')
%Generate the realization with the B coefficients
for j   = 1:N
	X   =  [X,Xt((p-1)*n+1:end)];
	Xn  = -A*Xt + iB0*randn(n,1);
	Xt  = circshift(Xt,n);
	Xt(1:n) = Xn;	
end


%Form the coefficient matrices
A = [eye(n) A];
B = [B0 B];
%free the symbol X
x = X';

%Generate the windowed estimate of the covariance matrix
[N,n] = size(x);
C = windowed_est(x, p);

%Build the sparsity pattern

clear Comp;
clear W;
clear burn;
clear X;
clear Xt;
clear Xn;
clear B0;
clear H;
clear stable;
clear me;
clear i;
clear j;

file = '../Data/cycle_graph_6_4_50.mat';
save(file)
fprintf('Results saved to %s\n',file)
%no need to do the p = k change, data is now lower case x with column per variable and row per sample.


